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Abstract 



In this paper, we extend the method proposed by Cochelin and Vergez in [A high 
order purely frequency based harmonic balance formulation for continuation of 
periodic solutions, Journal of Sound and Vibration, 324:243-262, 2009] to the 
case of non polynomial non-linearities. This extension allows for the computa- 
' ^ . tion of branches of periodic solutions of a broader class of nonlinear dynamical 

systems. 

■ The principle remains to transform the original ODE system into an ex- 

tended polynomial quadratic system for an easy application of the harmonic 
balance method (HBM). The transformation of non polynomial terms is based 
on the differentiation of state variables with respect to the time variable, shifting 
the nonlinear non polynomial non-linearity to a time-independent initial condi- 
tion equation, not concerned with the HBM. The continuation of the resulting 
algebraic system is here performed by the Asymptotic Numerical Method (high 
order Taylor series representation of the solution branch) using a further dif- 
ferentiation of the non polynomial algebraic equation with respect to the path 
parameter. 

A one d.o.f. vibro-impact system is used to illustrate how an exponential 
non-linearity is handled, showing that the method works at very high order, 
1000 in that case. Various kinds of nonlinear functions are also treated, and 
finally the nonlinear free pendulum is addressed, showing that very accurate 
periodic solutions can be computed with the proposed method. 
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1. Introduction 



In a previous issue, Cochelin et Vergez presented A high order purely frequency- 
based harmonic balance formulation for continuation of periodic solutions p|, 
targeting nonlinear dynamical systems. The technique addresses autonomous 
or forced nonlinear dynamical systems, using an arbitrary high order Harmonic 
Balance Method (HBM, see @, H|) together with the asymptotic numerical 
method (ANM, see [J, [B| for details) for the continuation. 

The proposed method automatically derives the HBM set of algebraic equa- 
tions to the desired order, provided that the governing equations are given in a 
quadratic formalism, ie, dynamical systems must be expressed as a set of first or- 
der nonlinear ordinary differential equations and nonlinear algebraic equations, 
all with purely quadratic polynomial nonlinearities at most (see eq. (3) in [![). 

The method has been used in various applications, such as the analysis of 
microbubbles in liquids by Pauzin et al. [6|, or the nonlinear oscillations of nano- 
and micro-electromechanical devices by Kacem et al. Q- A purely frequency- 
based stability analysis was also proposed by Lazarus and Thomas [8J, as a 
companion to the method presented in [ljj. 

Whereas classical harmonic balance technique is often limited to very few 
harmonics, the proposed method has the advantage of allowing an arbitrary 
high order of the Fourier series. The difficulty is the need to recast the orig- 
inal system of equations (which is generally not quadratic), into a quadratic 
framework without changing the nature of the problem, ie, the recast should be 
only a change of formlution not a change of the original problem. The recast 
of polynomial, square root, and rational functions has been treated in the cited 
paper, defining intermediate variables and using quadratic polynomial algebraic 
equations. The authors also set up the basis for the recast of a sine nonlinear- 
ity, as encountered in the simple, free pendulum equation 9"(t) + sm(8(t)) = 
(where the prime sign denotes time-derivation). They obtained a set composed 
of 4 nonlinear ordinary differential equations and 2 time-independent, nonlinear 
algebraic equations. But, the calculation was not carried out. 

In this paper, we propose a generalisation of the method proposed in [l[ 
for the treatment of ODE with non polynomial nonlinearities. This generalisa- 
tion greatly enlarges the field of application by overcoming the need for strictly 
quadratic reformulation. Therefore, the method can be used to compute peri- 
odic solutions families of most nonlinear dynamical systems. Notice that once 
a periodic solution is known, stability and bifurcation can be analyzed by using 
classical methods such as for instance, monitoring the eigenvalue of the mon- 
odromy matrix. This topic is not addressed hereafter. 

The present paper is organized as follows : in section [21 the case of the 
exponential function is discussed as an introductory example and the periodic 
solutions of a one d.o.f vibro-impact system are calculated using up to 1000 
harmonics (with MANLAB software) ; then, in section |3l the method is gen- 
eralized to a wider class of nonlinearities and in section |4j we show how this 
generalisation applies to the natural logarithm, trigonometric functions, as well 
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as non-integer power functions ; and last, in section [5j we use the proposed 
method to compute the periodic solutions of the simple free pendulum and we 
discuss the accuracy of the solution. 



2. An example using the exponential function: a regularised vibro- 
impact oscillator 

We consider a one-degree-of-freedom, mass-spring oscillator which is limited 
to the half-plane x < 1 by a rigid wall. The rigid wall reaction is modelled by 
an exponential function, with a coefficient a to tune the wall stiffness. 

Let x(t) denote the mass position, we look for periodic solutions of : 

x"(t) = -x(t) - lx'(t) - e^W" 1 ), (1) 

where I is a free parameter (see|2]for details on this equation). 

Before embarking in the treatment of ([TJ , we recall that two numerical meth- 
ods are used in [l[ : a high-order HBM technique and a high order Taylor series 
expansion for the continuation. The automation of these two techniques rely on 
the quadratic form of the nonlinearities. To comply with this framework, the 
method proposed by Cochelin et Vergez in [l[ uses the following formalism: 

m(Z'(t)) = c + l Cl + l (Z(t)) + lli(Z(t)) + q(Z(t), Z(t)) (2) 

where Z(i) is the unknown vector composed of time-dependent state variables, I 
is the continuation parameter, (co, Ci) are constant vectors, (lo, li, ni) are linear 
vector-valued operators, and q is a bilinear vector-valued operator. 

We will now present the successive recasts of system ([T]) needed to obtain the 
required general form ©. It should be noticed in passing that the quadratic 
form (0) has nothing to do with a second order Taylor expansion. The pas- 
sage form (fTJ) to (|2j) should not introduce any approximation, it is just another 
formulation of the same original problem. 

2.1. First-order recast 

We introduce an additional variable y{t) = x'(t) and rewrite equation ([T]) 

as: 

x'(t) =y{t) (3a) 
y'{t) = - x(t) - ly(t) - e^W')" 1 ) (3b) 

2.2. Quadratic recast of the exponential junction 

Introducing the additional variable e and its definition equation 

e(t) = exp a(x(t) - l) (4) 
we differentiate it with respect to the time variable to obtain 
e'{t) = ax' (t) exp a(x(t) — l) 

or equivalently 

e'(t) = ae(t)y(t) (5) 
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Note that for equation (O to be exactly equivalent to @, the following initial 
condition must be added: 



e(0) = exp a(x(0) - l) (6) 

We can now recast (15a| - (|3b|) into the following set of differential and algebraic 
equations: 

x'(t) =y(t) (7a) 
y'(t) = - x(t) - ly{t) - e(t) (7b) 
e'(t) =ae{t)y{t) (7c) 
e(0) = e a ^-^ (7d) 

At this point, the authors wish to underline the fact that the non polynomial 
nonlinearity has "moved" from an ODE equation (here: I3b[) to an algebraic 
equation (here: I7d|) which is not involved in the balance of harmonics: this is 
the key point of the method. 

2.3. Applying the harmonic balance method to the ODEs 

We will now apply the harmonic balance method (HBM) to the equations 

Cil-CTc]). 

Denote Z(t) = [x(t),y(t),e(t)] T , and the following operators: 

m(Z'(t))=[x'(t), y '(t),e'(t)] T 
c =ci =[0,0, 0] T 
lo(Z(t)) =[y(t),-x(t)-e(t),0] T 
h(Z(t)) = [0,-j/(i),0] T 
q(Z(t),Z(t)) =[0,0, ae(t) y (t)f 

The ODEs (|7all7c[) now take the required form: 

m(Z') = c + lCi + 1 (Z) + Ili(Z) + q(Z, Z). 



Eq. (|7all7bp are treated classically, using the method proposed by Cochelin 
and Vergez in by balancing for the H + 1 harmonics < h < H of Z(t), 
where H is the chosen truncation order of the Fourier series of the vector of 
(time-dependent) variables Z(i): 



H 

Z{t) = Z + (Zih-i cos{kujt) + Z 2h sm(kcjt)) (8) 

h=l 

For equation (|7c|) . only harmonics 1 < h < H need to be balanced. The 
reason is that the balance of the harmonic zero, which consists in equating the 
mean value of each side over a period, is not suitable here: both sides of the 
equation result from time-differentiation of a periodic quantity (namely e(t)), 
and has mean value which is therefore null. Thus, we only apply the HBM for 
harmonics 1 < h < H , where H is the truncation order. 
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The unknown (time-independent) variables are gathered in the state vector 

U: 

U = [{Zj}j=o..2H, A,w], 

whose size is N U =3(2H + 1) + 2. 

The algebraic equations resulting from HBM application to (|7all7bllTc]) , are 
then automatically put into the standard ANM quadratic formalism: 

L + L(U) + Q(U,U) =0 (9) 

with a constant vector Lo, a linear vectir- valued operator L, and a quadratic 
vector- valued operator Q. 

Now counting the number of algebraic equations: 



• 2(2H + 1) for the full HBM applied to f7a ) -([7b ]l 

• 2H for the partial HBM applied to (f7c|) 

• 1 for the time-independent equation (|7dl) 

• 1 for an additional phase equation (for instance: y(0) = 0) 

we reach a total number of N e =3(2H + 1) + 1 equations. 

Thus, we get N e — N u — 1, as is needed for a ID family of solutions. 

The reader shall notice that this algebraic system of N e equations for N u 
unknowns may now be solved by any continuation algorithm, such as for instance 
the very classical first order predictor with Newton corrector and pseudo arc- 
length parametrization. 

In the following, we aim at using an ANM continuation method because 
of its robustness as compared to predictor-corrector algorithm. However, the 
algebraic system is not fully quadratic, ie, equation @ is but not (|7d[) . hence, 
an additional step is necessary. 



2.4- Recast of the algebraic equation | Id) for ANM continuation 



We will now address the last recast that concerns non polynomial, nonlinear, 
algebraic equation (|7d[) . following the usual procedure presented in @ and 0. 
In the continuation process, the N u unknowns are function of a path parameter 
a. By differentiation of the non polynomial equation with respect to that path 
parameter a, quadratic equation can be obtained as shown below. 

For sake of clarity, let us define two new variables, eto and x t o, that represent 
e(0) and x(0) respectively: 

H H 

eto = eo + e 2 h-i x t p = x + }^x 2 h-i- 

h=l h=l 

where ei (respectively x\) denotes the i-th coefficient of the Fourier scries of e(i) 
(resp. x(t)) as defined for Z in the equation ©. 

By differentiation, equation (|7d[) is strictly equivalent to: 

Vo £ M, deto{a) = aeto{a) dxto(a) (10a) 
e*o(o=0) = e ^"^- 1 ) (10b) 
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where de — ^ denote differentiation with rexpect to a. 

Equation (|10a[) is now quadratic in (U, dU) and is suitable for an easy and 
efficient computation of the Taylor series used for the continuation. The reader 
is referred to appendix [C] for details on the formalism used to enter (|10aH10b"|) 
in the MANLAB software. 

2.5. Periodic solutions of the regularised vibro-impact 




Figure 1: Phase portrait of the regularised vibro-impact: family of periodic solutions 
in the phase plane (u,v). (a): non stiff regularisation (a = 20) with H = 50 harmonics, (b): 
stiff regularisation (a = 200) with H = 1000 harmonics. 

Figure [T] shows samples from the family of periodic solutions computed for 
two cases: a non stiff regularisation, with a = 20, using H = 50 harmonics (a) 
; and a stiff regularisation, with a = 200, using H = 1000 harmonics (b). In 
both cases, the continuation was run with a threshold of 10~ 10 on the residue. 

The phase portrait cycles are to be compared with those of a free, conserva- 
tive non-regular vibro-impact system, i.e. a wall modelled with an impact law 
using a restitution coefficient equal to unity : the family of periodic solutions is 
composed of origin-centered circles, when the amplitude is less than unity, and 
origin centered arc of circles closed by a vertical segment along the line x = 1, 
when the amplitude is higher than unity. 

3. General treatment of nonlinear functions 

Here, we discuss the general method for the recast of most nonlinearities 
into quadratic formulation. Note that the quadratic recast of rational functions 
has been given by Cochelin and Vergez in 

Let us consider a set of differential and algebraic equations: 

F(u,v,g(u)) = (11) 
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where v—u', F is at most quadratic in its arguments, and g is any nonlinear 
function of u. 



3.1. First order derivative 

We add a new variable w defined by w—g(u). By time-derivation, we obtain: 

w' = ^(u)v (12a) 
u;(0) = fl (u(0)) (12b) 

If dg/du can be written as a rational function of (it, v, w), then the equation 
x = dg/du can be recast into quadratic equations (possibly using additional 
variables). Consequently, the time-dependent equation (|12a|) becomes w' = xv, 
which is quadratic in x and v. 

We apply the HBM to this equation, but only for harmonics 1 < h < H, 
while the mean value (harmonic zero) will be constrained by the initial condition 
(I12b[) . The equation count is then equal to 2H +1, as in a standard HBM 
applied to a unique equation, which matches the number of variables : the 
2H + 1 coefficients of the Fourier series of w(t) up to harmonic H . In the case of 
autonomous systems, the period (or equivalently the angular frequency) is also 
unknown and one need to add a phase equation, as explained by Doedel in [l(| 
as well as Cochelin and Vergez in [if. 

If additional variables were to be used for the quadratic recast of x, all 
corresponding algebraic equations will be treated with full HBM, including the 
balance of the harmonic zero (the mean values). 

3.2. Second order derivative 

If x = dg/du cannot be expressed as a quadratic polynomial of the current 
variables, we differentiate it with respect to t: 

x' = 0(u)« (13a) 

*(0) = ^(u(0)) (13b) 

If d 2 g/du 2 can be expressed as a rational function of (u,v,w,x), then we 
define y=d 2 g/du 2 and the previous results apply: the time-dependent equation 
in (|13a[) is quadratic in y and v. 

We then have two differential equations, namely w'=xv and x'=yv, with 
two associated initial conditions. The differential equations are treated using 
HBM for harmonics 1 and higher, while the initial conditions will constrain the 
mean values of w and x (the nonpolynomial, nonlinear, algebraic equation is 
addressed by differentiation, as explained in section . 

If additional variables were used for the recast of the rational function y, all 
corresponding algebraic equations will be treated with full HBM, including the 
balance of the harmonic (the mean values). 

4. Recast of a few common non-polynomial nonlinearities 

For the quadratic recast of the exponential function, the reader is referred 
to section O 
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4-.1. Natural logarithm 

Given w denned as w(t) = \n(u(t)) (assuming u > 0). For sake of simplicity, 
we do not write explicitly the time dependence in the following. Differentiating 
with respect to the time variable, the definition equation becomes w' = u' /u, 
or equivalently, using x = w': 

( u' = xu 

\ w (t = 0) = ln(u(t = 0)) 
which is quadratic in u and x. 



4-2. Non-integer power 

Given w(t) = u{t) a where a £ I is a constant, then w' — au a ~ 1 u'. Using 
v = u' and x = w' , one gets: 

{ux — awv 
w(t = 0) = u(t = 0) a 

which is quadratic in (u,v,w,x). 



4-3. Trigonometric functions 

Given s(t) = sm(u(t)) and c{t) = cos(u(t)), we introduce v(t) — u'(t) and 
time-derivation of the definition equations of s and c gives: 

s' = cv 
c! = —sv 

s(t = 0) = sin(u(i = 0)) 
c(t = 0) = cos(u((i = 0)) 

which are obviously quadratic in (c, v) and (s, v) respectively. 

As for w(t) = tan(u(£)), time-differentiation leads to w' — (1 + w 2 )v. Using 
an additional variable x — 1 + w 2 , one gets the quadratic equation: 

w' = XV 

w(t = 0) = tan(w(i = 0)). 



5. Periodic solutions of the simple, free, nonlinear pendulum 

Denoting 8 the angle between the current position and the lower, rest posi- 
tion, the equation of motion of the free pendulum writes: 

6"{t) + sin [6(t)] =0 (14) 

Adding a damping parameter I, that will vanish along the family of periodic 
solutions, as explained in section [Al the equation of motion becomes: 

e"(t) + W'(t)+siii [6(t)] = (15) 
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Using three additional variables, v(i)=0'(t), s(t)=sin \0(t)] andc(i)=cos [0(i)l, 
we recast the equation (|T5|) into the following quadratic, differential-algebraic 
system: 



e'(t) = v(t) 


(16a) 


v'(t) = -s(t) - lv(t) 


(16b) 


s'(t) = c(t)v(t) 


(16c) 


c\t) = -s(t)v{t) 


(16d) 


s(0) = sin (61(0)) 


(16e) 


c(0) = cos (0(0)) 


(16f) 



Using the proposed method, we apply: 

• the full HBM to equations (I16aj|16bj) 

• the balance of the harmonics 1 and higher to equations (|16dll6d)) 

and since the frequency is also unknown, we add a phase equation: 

v(0) = 0. 



(17) 



Finally, the initial conditions (|16cM16fj) arc differentiated with respect to 
the path parameter. The equations are then in the right (quadratic) form for 
applying the ANM continuation, ie, computing a high orde Taylor series of the 
unknowns with respect to the path parameter: 

L +L(U)+Q(U,U) =0 
Lh(U) = Qh(U,U) 



The total number of equation is then: 
N e = 2(2Jf + l) + 2(2H 
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full HBM 



- ^_2_^ 

HBM /i>0 Phase cq. in Lh/Qh 



in LO/L/Q 

= 4(2# + l) + l 
while the number of variables is 

N u = 4(2JT+1) 



Fourier components of &,v,s,c 

4(2iJ+l) + 2 



Figure [2] shows the phase portrait and the amplitude- frequency diagram of 
the periodic solutions family of the nonlinear, free pendulum system, computed 
with H = 100 harmonics. The theoretical amplitude-frequency diagram of this 
system (dashes) corresponds to the following analytic formula: 



T(9 max ) = ^K(sin 2 (0 ma:c /2)), 



where K is the complete elliptic integral of the first kind (see [ll[ for instance). 
The computed curve (plaine line) is superimposed on the theoretical one. 
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Figure 2: Periodic solutions of the nonlinear free pendulum, (a): phase portrait (8,v) 
showing a few samples of the periodic solution family, (b): amplitude- frequency diagram, 
the computed curve is superimposed on that of the analytic formula. HBM using H = 100 
harmonics, ANM-threshold: 10~ 15 , the residue ||R.|| was kept under 10~ 14 . 



The branch was computed with 29 step of continuation. The computation 
was stopped when the relative error between the computed angular frequency 
and its theoretical value (given the amplitude) reached 0.1%, which occurs at 
the solution point: 

(0 ma x\ _ /0.999998tt rad\ 
\ u> J ~ 1^0.112801 rad/sj 

where the norm of the residual vector is | |R| |=2.76 10~ 15 . Increasing the number 
of harmonics only leads to a closer approach to the limit point (0 max ,oj)=(ir, 0), 
where the period is infinite^. 

Notice that two additional variables must be introduced for this simple one 
d.o.f examle. For a multiple d.o.f. system with many kinds of non-polynomial 
nonlinearities, the quadratic recast has to be applied to each individual non- 
linear term. This may require a great number of additional variables and the 
transformed system may be much larger than the original one. This is however 
the price to pay with this method, the more complex the original system, the 
larger the transformed quadratic system. 

6. Conclusion 

In this paper, we extended the works of Cochelin and Vergez presented in 
[H to the case of general nonlinearities. 

A vibro-impact with exponential regularization was presented, and its peri- 
odic solutions were computed with the proposed generalization. The case of a 
very stiff regularization demonstrated the capabilities of this numerical tool to 
deal with a very high number of harmonics in the harmonic balance method. 

Finally, we showed how to apply the method for the continuation of the 
periodic solutions of the simple, free, nonlinear pendulum. Our results confirm 
that the harmonic balance method with a high number of harmonics is both 
affordable and well suited to the ANM continuation framework. It has been 



1 This corresponds to a hctcroclinic connection between the two saddle- nodes 
(8 S N,Vsn) e {(-7T, 0), (tt, 0)} 



10 



successfully implemented in the MANLAB software which provide an interactive 
graphical user interface for the continuation, in a widely used programming 
environment. 

The principal limitation of the proposed method is due to the need of an 
algebraic relation between the nonlinear function and its derivative, or its prim- 
itives and the state variables. However, the numerous examples treated in this 
paper show a variety of such functions that might appear in nonlinear dynamical 
systems. 

Concerning future work, the companion frequency-based stability analysis 
presented in [8[ could be adapted to general nonlinearities, following the same 
approach. Also, automatic differentiation, such as the DIAMANT package pro- 
posed by Charpentier et al in [l2[ , could give a mean to simplify the user input 
to the bare nonlinear functions, letting the user free of any recast. 
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Appendix A Vibro-impact with exponential wall reaction 



A.l Model 

We consider a one-degree-of-freedom, mass-spring oscillator which is limited 
to the half-plane x < 1 by a rigid wall, where x(t) denotes the position of the 
mass. 

The rigid wall reaction is modelled by an exponential function, with a coef- 
ficient a to tune the wall stiffness : 

F w {x) = -e a{ - x ~ l \ 

Thus, the regularised vibro-impact system is governed by the following equa- 
tion: 

x"(t) = ~x(t) - e"^*)- 1 ), (A.l) 

where the prime sign denotes time-differentiation. The force F w (x) that reflects 
the wall effect derives from a potential energy so that problem (|A.1[) keeps the 
property of being energy-conservative. 

A. 2 Dissipative recast for continuation 

Muhoz-Almaraz et al. [IH showed that, in conservative Hamiltonian sys- 
tems, periodic orbits generally belong to a one-dimensional family of periodic 
solutions, parametrised by the value of the first integral (here, the total energy), 
which is not an explicit parameter of the system. To compute this family of 
periodic solutions in the standard continuation framework R(x{t),l) = 0, we 
perturb the initial equation with a damping term added to the right-hand side 
of (|A.1[) . The system is then embedded into a general, dissipative system: 

x"{t) = -x(t) - lx'{t) - e^W" 1 ), (A.2) 

where I is a free parameter of the continuation. 

The perturbed system (|A.2|1 possesses periodic solutions that are exactly 
those of the unperturbed, conservative system (lA.ip . if and only if I = 0. This 
way, the additional parameter I allows us to compute the periodic solutions 
of the conservative system (|A. 1 1) using the classical framework for dissipative 
systems possessing an explicit control parameter. 

Appendix B Classical ANM: quadratic framework 

The reader is referred to Cochelin and Vergez [l[ (sections 2.4 and 2.5, 
pp. 248-250), for the details concerning the principle and the implementation 
of the ANM in the classical, quadratic framework. 

Appendix C Extended ANM framework 

C. 1 Principle of the series computation 

Given a nonlinear system f(U) = with Nd equations and N u unknowns, 
whose differentiated form reads: 

Lh(TJ) = Qh(U,U) (C.l) 
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where Lh : -> R Nd is a linear, vector-valued operator ; Qh : R w " x K w " ->■ 
Br*" 1 is a bilinear, vector- valued operator. 

Assuming a known regular solution Uo of this system, we write the branch 
of solutions passing through this point as a Taylor series expansion: 

U(o) = U + 0U1 + a 2 U 2 + a 3 U 3 + • • • + a"U n . (C.2) 

where the branch is parametrised using the pseudo-arclength parameter a de- 
fined as: 

o = CU-U )*Ui. (C.3) 

Differentiating U reads: 

U = {Ui + 2aU 2 + 3a 2 U 3 + ... + raa n-1 U n }da. (C.4) 

Then, substituting both (|C.2I) and (IC.4I) in system (IC.1[) and equating each 
power of a (up to order n) to zero gives: 

• power 0: Lh(Ui) — Qh(U , Ui) = 0, which can also be written Jhu -Ui = 
where Jhu € K ArdXArd+1 is the jacobian matrix of / evaluated at Uo- This 
linear equation in Ui thus gives the term of order 1 of (|C.2p . 

• power 1 < p < n - 1: Jh Uo .U p+1 = £f =1 £±^pQh(Ui, Up+i_<). This 
linear equation gives the term of order p + 1 of (|C.2I) . 

The original nonlinear problem is thus replaced by a cascade of n linear systems, 
which all share the same matrix: Jhu ■ 

However, at each order, the linear systems have Nd + 1 unknowns and only 
Nd equations. For each linear system, the additional equation is obtained by 
substituting the series (|C.2I) into the definition (|C3[) of the path parameter a: 

• order 1 : UjUi = 1 

• order 2 < p < n : XJ\U p = 

C.2 Implementation in MANLAB: the example of the vibro-impact 

As for the classical framework, the only user input to the MANLAB software 
consists in M-functions for the operators Lh and Qh, as well as for / (for the 
residue computation only), and a starting point Uo. 

In the case of the vibro-impact system presented in section[21 some equations 
are quadratic (those resulting from the HBM and those for the definition of e f o 
and xto) while the last one is not. We thus separate the equations in two parts: 
those resulting from the HBM, that will appear in the L0, L and Q operators, 
and the last one, that will appear in the Lh, Qh operators. 

For the present problem, the state vector is: 

U = ( xq, 2/q, eg , xi,yi,ei , x 2 , yi,e 2 , 

Zo Zi Z2 

X2H-l,y2H-l,e 2 H-l,X 2 H, V2H , e 2 H,l, W, &tO, X t o) 

Z2H-I Z2if 

and its size is N u = 3(2H + 1) + 4. 
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For the quadratic part of the problem, the subsystem contains the following 
number of equations: 

N e = 2(2H + 1) + ^H_, + 2 + = 3(2# + 1) + 2. 

full HBM HBM h^O dof. of e t0 ,x ta phase eq. 

The content of LO . m, L . m and Q . m are not listed here, as it is the direct result of 
the harmonic balance as presented in [l[ . These three functions return a vector 
of size N e — 1. 

For the non quadractic part, the subsystem contains iVj = 1 equation. The 
content of Lh.m, Qh.m and f .m is listed below, with a = 200: 

function [Lh] = Lh(dU) 
Lh=zeros(l,l) ; 
Lh = dU(end-l) ; 



function [Qh] = Qh(U,dU) 

Qh = zeros(l,l) ; 

Qh = 200*U(end-l)*dU(end) ; 



function [f] = f (U) 

f = zeros (1 , 1) ; 

f = exp(200*(U(end)-l)) ; 
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List of Figures 



1 Phase portrait of the regularised vibro-impact: family of 
periodic solutions in the phase plane (u, v). (a): non stiff regular- 
isation (a = 20) with H = 50 harmonics, (b): stiff regularisation 
(a = 200) with H = 1000 harmonics 

2 Periodic solutions of the nonlinear free pendulum, (a): 
phase portrait v) showing a few samples of the periodic so- 
lution family, (b): amplitude- frequency diagram, the computed 
curve is superimposed on that of the analytic formula. HBM us- 
ing H = 100 harmonics, ANM-threshold: 10~ 15 , the residue ||R|| 
was kept under 10 -14 
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